/* http://shootout.alioth.debian.org/u32q/benchmark.php?test=fasta&lang=gcc&id=4 */
/* adapted from C GNU gcc sample */

#define WIDTH		60
#define MIN(a,b)	((a) <= (b) ? (a) : (b))
#define NELEMENTS(x)	sizeof(x)

auto n = 1000;

struct aminoacid {
    float32_t	p;
    int8_t	c;
};

auto IM = 139968;
auto IA = 3877;
auto IC = 29573;
auto last = 42;
auto myrandom(max) {
    last = (last * IA + IC) % IM;
    return max * last / IM;
}

void accumulate_probabilities(aminoacid genelist[], len) {
    auto cp = 0.0;
    auto i;
    for (i = 0; i < len; ++i) {
	cp += genelist[i].p;
	genelist[i].p = cp;
    }
}

/* This function prints the characters of the string s. When it */
/* reaches the end of the string, it goes back to the beginning */
/* It stops when the total number of characters printed is count. */
/* Between each WIDTH consecutive characters it prints a newline */
/* This function assumes that WIDTH <= strlen (s) + 1 */
void repeat_fasta(s[], count) {
    uint8_t buf[] = new uint8_t[WIDTH + 1];
    auto off = 0;
    do {
	auto line = MIN(WIDTH, count);
	auto pos = 0;
	do {
	    if (off >= sizeof(s))
		off = 0;
	    buf[pos++] = s[off++];
	} while (pos < line);
	buf[line] = '\n';
	print("%*s", line + 1, buf);
	count -= line;
    } while (count);
}

/* This function takes a pointer to the first element of an array */
/* Each element of the array is a struct with a character and */
/* a float number p between 0 and 1. */
/* The function generates a random float number r and */
/* finds the first array element such that p >= r. */
/* This is a weighted random selection. */
/* The function then prints the character of the array element. */
/* This is done count times. */
/* Between each WIDTH consecutive characters, the function prints a newline */
void random_fasta(aminoacid genelist[], count) {
    uint8_t buf[] = new uint8_t[WIDTH + 1];
    do {
	auto line = MIN(WIDTH, count);
	auto pos = 0;
	do {
	    auto r = myrandom(1.0);
	    auto i = 0;
	    while (genelist[i].p < r)
		++i; /* Linear search */
	    buf[pos++] = genelist[i].c;
	} while (pos < line);
	buf[line] = '\n';
	print("%*s", line + 1, buf);
	count -= line;
    } while (count);
}

aminoacid iub[] = {
    { 0.27, 'a' },
    { 0.12, 'c' },
    { 0.12, 'g' },
    { 0.27, 't' },
    { 0.02, 'B' },
    { 0.02, 'D' },
    { 0.02, 'H' },
    { 0.02, 'K' },
    { 0.02, 'M' },
    { 0.02, 'N' },
    { 0.02, 'R' },
    { 0.02, 'S' },
    { 0.02, 'V' },
    { 0.02, 'W' },
    { 0.02, 'Y' }};

aminoacid homosapiens[] = {
    { 0.3029549426680, 'a' },
    { 0.1979883004921, 'c' },
    { 0.1975473066391, 'g' },
    { 0.3015094502008, 't' }};

accumulate_probabilities(iub, NELEMENTS(iub));
accumulate_probabilities(homosapiens, NELEMENTS(homosapiens));

uint8_t alu[] ="\
GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG\
GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA\
CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT\
ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA\
GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG\
AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC\
AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA";

print(">ONE Homo sapiens alu\n");
repeat_fasta(alu, 2 * n);
print(">TWO IUB ambiguity codes\n");
random_fasta(iub, 3 * n);
print(">THREE Homo sapiens frequency\n");
random_fasta(homosapiens, 5 * n);
